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1 Documentation 

Hydra is an adaptive particle-particle, particle-mesh plus smoothed particle hydro- 
dynamics N-body simulation program. It can be used with either periodic or isolated 
Q\ ■ boundary conditions. A compiler flag allows the gas calculation to be turned off, 

! converting Hydra to a collisionless mode. 

Q-f An installation guide comes with the release which provides information about 

how to set up and compile Hydra with parameters appropriate for your particular 
problem. This documentation describes the Hydra input files in more detail. 

In what follows code extracts and variables will be written in typewriter font. 
Program and subroutine names will appear in bold face and directory names will be 
bracketed <thus>. 



1.1 Data format 

The structure of data files is defined in dumpdata.f: 

real rm(N) ,r(3,N) ,v(3,N) ,dn(N) ,e(N) ,h(N) 

integer itype(N) 

write (8) ibuf , ibuf 1 , ibuf 2 

write (8) rm 

write (8) r 

write(8) v 
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write (8) h 
write (8) e 
write (8) itype 
write (8) dn 

(i) The first line writes out the values of important parameters and variables as 
listed in pinfo.inc: 

common/par am/ it ime , it stop , itdump , itout , time , atime ,ht ime , dtime , 
& Est , T , Th , U , Radiat i on , Esum , Rsum , cpu , 

& t stop , tout , icdump , padding , Tlost , Qlost , Ulost 

common/start/ irun , nob j , ngas , ndark , L , intl , nlmx , perr , 
& dtnorm, sf tO , sf tmin , sf tmax ,hl00 ,boxl00 , zmet , spcO , 

& lcool , rmgas , rmdark , rmnorm , t start , omegaO , xlambdaO , hOtO 

common/outputtime/toutl 

real toutl(25) 

dimension ibuf 1(100) ,ibuf2(100) ,ibuf (200) 
equivalence (ibuf 1 , itime) , (ibuf 2, irun) , (ibuf ,toutl) 

Each of these variables serves the following function; 
itime, itstop, itdump, itout: step number, stopping step, dumping interval for 
backups, not used 

time, atime, htime, dtime: time, expansion factor, hubble parameter, timestep 
(all in code units) 

Est, T, Th, U: starting, kinetic, thermal and potential energies in code units 
Radiation, Esum, Rsum, cpu: radiated energy, energy integral, total radiated en- 
ergy, cputime counter 

tstop, tout, icdump, padding: not used, next dump time, position in list of 
dump times, padding for isolated boundary conditions 

Tlost, Qlost, Ulost: lost kinetic, thermal and potential energies in grid units 
(isolated boundary conditions only) 

irun , nob j , ngas , ndark : run number, total number of particles, number of gas, 
dark particles 

L, intl, nlmx, perr: grid size (top level), interlacing switch, max number of re- 
finement levels, maximum 2-body percentage force error 

dtnorm, sftO, sftmin, sftmax: timestep multiplier, current-day softening (grid 
units), min, max softening (grid units) 

hlOO, boxlOO, zmet, sftO: hubble parameter, boxsize (/i _1 Mpc), metallicity of 
gas, average particle spacing 

lcool, rmgas, rmdark, rmnorm: cooling switch, mass of gas, dark matter particle, 
force normalisation 
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t st art, omegaO, xlambdaO, hOtO: start time, omega , lambda , H t (last three 
for expanding box only) 
toutl: list of desired output times, 
(ii) data arrays 

Note that the following arrays are defined for all particles even though some of 
them are irrelevant for dark matter particles. This is to simplify book-keeping. We 
recommend setting e, h and dn to zero for dark matter particles. 

Array used for 

rm mass 

r positions 

v velocities 

e thermal energy 

h sph smoothing length 

dn density 

itype particle type (meaning set in itype.inc) 
The current settings for itype are: 



itype 


used for 


-2 


particle does not exist (useful for merging particles) 


-1 


star (equivalent to dark matter) 





dark matter 


1 


gas 


2 


temporary setting for an already completed gas particle 



sph evaluation 

These can be changed to label particles in different ways but you should check 
that Hydra handles the resultant types correctly. 

1.2 Units 

WARNING: the units used internally by Hydra differ from those used in the datafiles. 
The reason for this is that internally Hydra uses positions in the interval [1,L+1), 
where L is the number of FFT grid cells across the box. The data is either mapped 
into this interval exactly or, for isolated boxes and refinements, a padding region is 
left around the outside of the box. L can be different in different refinements and 
can vary from run to run - it has no physical meaning but is purely a computational 
device. The units also differ in isolated and periodic simulations: 
(i) isolated simulations (usually compiled with -DISOLATED). 

• length: positions run from to 1 in each co-ordinate direction. The length-unit 
is boxlOO h^Mpc. 

• mass: currently set to 10 10 M Q . Can be altered. 



only meaningful for gas 
only meaningful for gas 
only meaningful for gas 
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• time: currently set to 10 10 yr. Can be altered. 

• speed: length/time 

• density: number of particles per box volume. Converted internally to a number 
density of ions + electrons using nunit. 

• temperature: stored in units of internal energy (ie speed 2 ). Convert to ergs 
by multiplying by eunit = vunit 2 ; convert to temperature by multiplying by 
Kunit = eunit*2 . *mum/3 . /kb 

(ii) cosmological simulations (periodic boxes) 

• length: positions run from to 1 in each co-ordinate direction. The length-unit 
is atime*boxlOO h^Mpc. 

• mass: the total mass in the box is normalised to match the desired value of Q,q. 
Thus particle masses are all relative. 

• time: normalised to unity at the present, ie tunit=l at the end. t is calculated 
automatically from input cosmological parameters. 

• speed: The output velocities are dx/dt where x is the spatial coordinate in the 
range [0,1). The peculiar velocity is thus atime * dx/dt and the proper velocity 
datime/dt*x + atime*dx/dt. 

• density: number of particles per box volume. Converted internally to a number 
density of ions + electrons using nunit. 

• temperature: stored in units of internal energy (ie speed 2 ). Convert to ergs by 
multiplying by atime 2 *eunit = vunit 2 ; convert to temperature by multiplying 
by Kunit = atime 2 *eunit*2 . *mum/3 . /kb 

The units are initialised in the routine inunit.F to which you should refer if in 
doubt. They are also written at the head of the log file. 

The magnitude of the force-scaling depends upon the system of units which is in 
use. It is set by the factor rmnorm defined at the head of updaterv.F. 



1.3 Creating initial conditions 

The input data arrays which need to be defined for all particle species are rm, r, v 
and itype. In addition gas particles need an initial temperature, e, and sph smooth- 
ing length, h. This latter quantity may be estimated as the mean interparticle sepa- 
ration, dn -1 '/ 3 '). dn itself is not required. 

For dark matter runs the following are the minimum set of input parameters which 
need to be defined: itime, time,irun, nobj, hlOO, boxlOO, omegaO (periodic boxes 
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only), xlambdaO (ditto), toutl. In addition for gas runs the cooling flag needs to be 
set: lcool=l for cooling, otherwise, and for lcool=l then the metallicity, zmet, 
should be defined in terms of the Solar value. 

An example initial conditions generator for an isolated box is createtop.f, included 
with this distribution. A periodic initial conditions generator is also included, the 
cosmic package. 



1.4 Running Hydra 

Hydra works within and below a runtime directory. An example of the required 
directory setup was automatically built by unpacking the tar file. The executable 
(hydra) will be copied to the run directory <rundir> automatically when it is made 
(the destination directory can be reset by changing RUNDIR in the makefile). When 
Hydra is executed the runlog (hydra > runlog) and a summary file (pr????.log 
where ???? is the run number) are also written out here whilst the data is written to 
(and read from) a sub-directory <data>. The Green's functions the code uses to do 
its PM calculation are saved to the sub-directory <greenfn>. If you have more than 
one <rundir> then it is usually best to link together the greenfn directories by using 
In -s .... So the directory tree looks like this; 

<rundir> < hydra, prun.dat, logfiles 

/ \ 
/ \ 
/ \ 

datafiles — > <data> <greenfn> < — Green's functions 



Hydra uses a short file to set its basic run parameters. This file, <rundir>/prun.dat 
has the following format: 

d3712.0857 startup data file name 

3724 irun; run number for this run 

858 10 itstop itdump 

1.0 dtnorm 

0.02 sft(Mpc/h) 

1 refinements on/ off 

The first line contains the name of the datafile you wish to start with. This is built 
up of d????.nnnn where ???? is a run number and nnnn is the step number. You 
can reset the run number on the second line. Here we are starting with the 857th step 
of run 3712, calling this run 3724. The third line contains the step number you wish 
to stop at and how frequently you require incremental backup files, dtnorm is the 
timestep normalisation: this should normally be set to unity (see Section 3.8 below). 
Next comes the Plummer softening in units of h~ 1 Mpc (see Section 3.7 below). The 
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final line contains a switch for turning refinement placing on or off. You can prevent 
Hydra placing refinements and convert it into a P3M code by setting this parameter 
to 0. 



1.5 Log files 



Hydra sends output to the screen and to a summary file. The screen output can be 
redirected to a log file (e.g. nohup hydra > run????. log ) but it is quite large and 
useful mainly for diagnostic purposes. 

The summary file is called pr????.log where ???? is the run number. One line is 
written to the summary file every step. Its form differs slightly between periodic and 
isolated runs: 

(i) periodic boxes 



step cputime time 

857 927.48 0.5440 

858 966.93 0.5445 



redshift 
0.5005 
0.4997 
W 

7.41E+07 
7.42E+07 



K 

5.00E+07 
5.01E+07 
error 
0.0065 
0.0065 



U 

2.52E+06 
2.52E+06 
(K + U)/a 
7.89E+07 
7.88E+07 



W/a 
1.11E+08 
1.11E+08 



where K— Kinetic energy, U— Thermal energy, W= Potential energy and a— 
atime. The final two columns should be in the ratio 2:3 (approximately - ignoring 
softening) and should remain constant in time during the linear regime. The error is 
defined in terms of the ratio of the change in the Layzer-Irvine energy integral, /, to 
the potential energy. / is defined as I = K + U + W + J [2 (K + U)+ W]da/a + J Ldt 
where L=radiated power and t = time. See Couchman, Thomas & Pearce (1995) for 
details. 

(ii) isolated boxes 



step cputime 

112.58 

1 112.25 



time 
0.0000 
0.0438 



K 

0.00E+00 
6.12E-03 
error 
0.0000 



U 

0.00E+00 
3.20E-10 

lost K 
0.00E+00 



W 

1.26E+01 
1.26E+01 

lost U 
0.00E+00 



lost W 
0.00E+00 



0.0005 0.00E+00 0.00E+00 0.00E+00 

The error estimate in this case is [K + U + W — lost(K + U + W) — starting (K + 
U + W)]/U. 

The screen output begins with some warning messages informing the user if any 
parameters in the input data file have been altered. It also prints the size, L, of the 
top-level grid and confirms the presence or absence of refinements. Next follows a list 
of units (remember to include appropriate atime factors in expanding boxes). Next 
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comes some diagnostic information: 



itime, time, atime, htime: 1295 1. 


000 I." 


r>AA 

JUU 


U . 00 I 






refinement: 0, scaling factor: 1. 


00 










ndark, ngas, ndone, nstar: 32768 32768 




A 

u 


U 




noveredge, novergrid: 2578 


3697 










gravity: av,min,max # neighbours 


110.22 




u 


loo4 




sph: av,min,max # neighbours 


23.48 




u 


101 




av,min,max density (ref.) 7. 


6971 


0, 


.0139 


383. 


.5992 


av,min,max energy (base) 50. 


3896 


0, 


.0000 


45135. 


.1523 


refinement: 1, scaling factor: 3. 


29 










ndark, ngas, ndone, nstar: 2247 906 




1241 







noveredge, novergrid: 213 


94 










gravity: av,min,max # neighbours 


41.12 







391 




sph: av,min,max # neighbours 


43.19 




28 


96 




av,min,max density (ref.) 1. 


0347 


0. 


.0239 


7. 


.8099 


av,min,max energy (base) 124. 


5373 





.0739 


433. 


.4357 



...similarly for other refinements... 
timestep = 2 . 6515502E-04 (2 . 6515502E-04 7 . 3740509E-04 9 .3750297E-02) 

(a) 'scaling factor' is the expansion factor of the refinement grid relative to the 
base level. 

(b) ndark, ngas, ndone and nstar give the number of particles of each type within 
that refinement: ndone refers to those gas particles whose sph forces have already 
been calculated at the previous level. 

(c) The average minimum and maximum number of gravity neighbours refers to the 
pp force calculation: the average value should be about 100 for efficient distribution 
of work between PM and PP. 

(d) The minimum number of sph neighbours can drop below 32 in low-density 
regions because we only search for neighbours out to a radius of approximately 2.2 
times the grid-spacing: within refinements the minimum number of sph neighbours 
should be close to 32. The maximum number of sph neighbours can be larger because 
the smoothing length can never drop below the gravitational softening. 

(e) The average, minimum and maximum density of the gas particles is given in 
refinement units (ie particles per refinement grid cell). The average, minimum and 
maximum temperature of gas particles is scaled back to the base level units. These 
units both differ from those in the data files. These entries are omitted from the log 
if there are no gas particles within that refinement. 

(f) The final line gives the timestep and the constraint on the timestep from 
accelerations, velocities and the hubble expansion (see Section 3.8 below). 
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1.6 Softening 

Hydra uses softening shape function which has compact support: 

SUBROUTINE S2(r,a,grav,dpot) 
xi=2.*r/a 

if(xi.ge.O. .and. xi . It . 1 . )then 

grav=xi* (224 . +xi*xi* (-224 . +xi* (70 . +xi* (48 . -xi*21 . ) ) ) ) /35 . /a**2 

dpot=(208.+xi*xi*(-112.+xi*xi*(56.+xi*(-14. 
& +xi*(-8.+xi*3.)))))/70./a 
else if(xi.ge.l. .and. xi . It . 2 . )then 

grav= (12 . /xi**2-224 . +xi* (896 . +xi* (-840 . +xi* (224 . +xi* (70 . + 
& xi*(-48.+xi*7. ))))))/35./a**2 

dpot= (12 . /xi+128 . +xi* (224 . +xi* (-448 . +xi* (280 . 
& +xi* (-56 . +xi* (-14 . +xi* (8 . -xi) )))))) /70 . /a 

else 

grav=l . /r**2 

dpot=r*grav 
end if 
RETURN 
end 

It is approximately equivalent to a Plummer softening with extent a/2.34 (the 
central value of r/force is the same in each case). The user supplied value for the 
softening in prun.dat is assumed to be a Plummer softening and is therefore multiplied 
by 2.34 before converting to grid units. If the ISOLATED or COMOVING flag is set 
during compilation then the softening is held fixed. Otherwise it scales with time as 
soft = min(0.6, sftO/atime) where sftO is the user-supplied value and the maximum 
value of 0.6 is a numerical constraint. 

When choosing a softening it is important to strike a balance between the compet- 
ing desires of high spatial resolution and long timesteps. dtime scales approximately 
in proportion to t s = \J (soft 3 /Gm) where m is the mass of the highest-mass objects 
in the simulation. One should also be aware of the danger of artificial two-body re- 
laxation which will occur in high-density regions on a timescale t 2 = V A5Nt s where 
iV is the number of particles within the softening. See Thomas & Couchman (1992) 
for more more a detailed discussion. 

1.7 Timestepping 

Our choice of timestepping algorithm is described in Couchman, Thomas & Pearce 
(1995). We evaluated several algorithms and settled for a simple PEC scheme which 
utilises only the current positions (plus velocities and temperatures for gas particles) 
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to determine the forces. This has the advantage of keeping storage to a minimum 
and also allowing arbitrary changes in timestep if the force changes abruptly in time 
(as it can do in the vicinity of shocks). The scheme is equivalent to Leapfrog (but 
less memory-efficient) if there are no gas forces. 

The timestep used by Hydra is not spatially variable but it does vary from one step 
to the next. The timestep is set by examining the greatest acceleration and velocity 
present at the current time: 

dta=(sof t2/asqmax) **0 . 25 
dtv=sqrt (sof t2/vsqmax) 
dtime=min(0.25*dta,0.4*dtv) 
if (htime.gt.le-30) then 

dtime=min(dtime , . 0625/htime) *dtnorm 
else 

dtime=min(dtime , 1 . ) *dtnorm 
end if 

where sof t2 is the softening length squared, asqmax is the square of the maximum 
acceleration and vsqmax is the square of the maximum velocity. For expanding boxes 
the htime condition ensures that cooling due to the Hubble expansion is followed 
accurately: this dominates at early times. 

1.8 SPH parameters 

The sph smoothing length is chosen so as to encompass approximately 32 neighbours. 
Some people prefer a higher value but we have found that this makes an impercep- 
tible difference to the results (32 particles gives significant shot-noise scatter in the 
interpolated density for a Poisson distribution but the actual particle distribution is 
much more uniform than this). In low-density regions the maximum sph search length 
(which is limited to approximately 2.2 times the grid-spacing) encloses fewer neigh- 
bours; in high-density regions the minimum sph search-length (which is set equal to 
the S2 softening parameter — see Section 3.7 above) may enclose more particles. If a 
particle's sph search length either extends beyond the refinement boundary or is more 
than 2.2 grid-spacings in size then the sph force for that particle is calculated at the 
previous refinement level (resulting in a loss of efficiency). The number of particles 
for which this is done is given in the log. 
We use a compact smoothing kernel: 

wnorm=l . /4. /pi 

if (x.le.1.0) then 

kernel=wnorm* (4.-6. *x**2+3 . *x**3) 
else if (x.le.2.0) then 
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kernel=wnorm* (2 . -x) **3 
end if 

where x=r/h is the ratio of the particle separation to the smoothing length, except 
that the gradient is modified so as to ensure a force which decreases monotonically 
with radius (see kernel. f and Thomas & Couchman 1992). 

The artificial viscosity is based on the bulk divergence of the flow, not on the 
relative pairwise velocity of particles. Schematically: 

f rc=twothird*e 

if (divpph. It . . ) f rc=f rc+(bvisc*divpph-avisc*cs) *divpph 

where divpph=h*div(v), cs is the sound speed and the parameters avisc and bvisc 
are set to 1 and 2, respectively. The actual code, which includes a correction for the 
hubble expansion, is buried in the heart of shgravsph.F. 

1.9 Cooling 

Cooling can be turned on or off by setting the parameter lcool (0=off, l=on). The 
supplied cooling function is a simple series of power-law fits to the optically-thin 
radiative cooling code of Raymond, Cox & Smith: 

if (t.lt.le5) then 

Lambda= ( 1 . 4e-28+ 1 . 7e-27*zmet ) *t 
else 

Lambda=5 . 2e-28*t**0 . 5+1 . 4e-18*t** (-1 . ) +1 . 7e-18*t** (-0 . 8) *zmet 
end if 

where t is the temperature in Kelvin. This fit contains contributions from both 
bremsstrahlung and recombination cooling from hydrogen and helium plus a variable 
contribution from metal lines with an assumed abundance of zmet times solar. The 
fit is good to within about a factor of two above lO 4 ^ - - below this temperature 
the cooling function drops precipitously and we take it to be zero. We intend to 
incorporate a more accurate cooling function in a future release. 

Because the cooling time is often shorter than the dynamical time, we do not use 
it to limit the timestep. Instead we allow the particles to evolve adiabatically, then 
cool them, assuming a constant density, at the end of the timestep. 

1.10 Brief description of the code 
MAIN: 

startup: read in parameters and data 
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inunit : define the units 
loop until finished: 
updaterv: PEC step 

I accel: acceleration including hubble drag; timestep evaluation 
I force: acceleration evaluation 

I rfinit: initialise refinements 

I ref force: see below 

I clist: create particle lists for refinements 

I loop over refinements: 

I load: load particles into refinement 

I | ref force: see below 

I | uload: unload particles from refinements 

I end loop 

I infout: write summary file 

I output: write out data and backup files 

end loop 

end 



ref force : 



green : evaluate, or read, green's function 

mesh : evaluate PM accelerations 

list : sort particles into search cells 

refine : determine the position of subref inements 

shforce : tabulate PP force and potential 

shgravsph : apply PP and SPH forces; 

write out diagnostic information 



1.11 Cosmological initial conditions 

A cosmological initial conditions generator - cosmic comes with the Hydra pack- 
age. It produces an initial conditions file in Hydra format from a supplied data file, 
cosmic.dat. There are in act 3 executables in the cosmic package. These are cos- 
mic itself which takes cosmic.dat for input and outputs the Hydra datafile into 
<data>d????.0000, where ???? is irun from cosmic.dat and r.pert is a pertur- 
bation file read by peakfindfft. peakfindfFt searches r.pert for peaks of a given 
size (the default is 8/i _1 Mpc) and outputs the position and 5 of the biggest peak into 
peak.dat. Finally there is pnsum which compares the sums of waves in boxes to 
the true power spectrum in spheres of radius 8h~ 1 Mpc (only useful if boxsize exceeds 
this) using the parameter file cosmic.dat. 
cosmic.dat has the following format; 
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1000 
4096 



-1066601 



irun, iseedl 
ndark 



20 
1. 
.06 
2 



0.5 
0. 
0. 




boxl00,hl00 
omegaO , xlambdaO 
omegabO , gamma 
ispec, ind 
atime , sigma8 
lcool,zmet 



0.02 




0.6 
0.5 



Notes; 



• ndark must be the cube of a power-of-2 

• omegab=0 for dark matter only, omegab> gives ndark each of dm and gas with 
correct mass ratio 

• if gamma<= then uses formula in Viana & Liddle (1996) 

• — ipsec = 1 - power law (index ind) 

— ipsec = 2 - cdm 

— ipsec = 3 - hdm 

• sigma8 is normalised to initial time using formulae in Viana & Liddle (1996) 

• lcool=l for cooling (0 for no cooling), zmet is in solar units 

1.12 Utility programs 

1.12.1 readheader 

The program readheader tells you some useful information about your initial con- 
ditions (or any dataset). It will automatically be placed into <rundir> when made, 
readheader prompts you for the name of a dataflle and checks to see if you have set 
Nmax and Lmax correctly and also shows the ranges which your data spans if you 
request further information. This utility is very useful for checking data consistency. 

1.12.2 hydra2tipsy 

hydra2tipsy converts Hydra output files to tipsy format. It converts positions to 
h^Mpc, velocities to km/s, temperatures to Kelvin and densities to overdensities 
(or atoms per cm 3 if an isolated box has been asked for). hydra2tipsy is also 
automatically placed into <rundir> by make and prompts for a Hydra datafile. The 
output file is in the format tip????.xxxx where ???? is the run number and xxxx 
is the current step. Tipsy, an N-body visualisation package can be obtained from; 

• http://www-hpcc.astro.washington.edu/tools/TIPSY/ 
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